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Abstract. We study a driven-dissipative array of coupled nonlinear optical 
' (— I '. resonators by numerically solving the Von Neumann equation for the density matrix. 

[ We demonstrate that quantum correlated states of many photons can be generated 

also in the limit where the nonlinearity is much smaller than the losses, contrarily 
to common expectations. Quantum correlations in this case arise from interference 
between different pathways that the system can follow in the Hilbert space to reach its 
steady state under the effect of coherent driving fields. We characterize in particular 
two systems: a linear chain of three coupled cavities and an array of eight coupled 
cavities. We demonstrate the existence of a parameter range where the system emits 
photons with continuous- variable bipartite and quadripartite entanglement, in the case 
of the first and the second system respectively. This entanglement is shown to survive 
realistic rates of pure dephasing and opens a new perspective for the realization of 
quantum simulators or entangled photon sources without the challenging requirement 
of strong optical nonlinearities. 



PACS numbers: 42.50.-p, 71.36,+c, 42.50.Ex 
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1. Introduction 

The physics of interacting Bose particles on a lattice offers a general framework within 
which several systems in condensed matter, atomic and optical physics are modelled. 
In the field of ultracold atoms in particular, the striking experimental demonstrations 
of phases with strong quantum correlations - such as the Mott insulator phase [H El |3] 
- and of the associated phase transitions, have convincingly opened the way to the 
physical realization of quantum simulators [H [5] . A similar perspective [6] has vigorously 
stimulated the fields of quantum optics and cavity quantum electrodynamics (CQED), 
supported by the technological advances in the fabrication of optical micro- and nano- 
resonators with exceptionally long photon lifetimes. In this context, a major objective 
consists in the implementation of the physics of strongly correlated bosons in a system 
of photons confined in an array of coupled optical resonators. In such a system, the 
effective on-site interaction between photons can originate from optical nonlinearities, 
due to some nonlinear medium embedded in each cavity. In particular, two typical kinds 
of nonlinearities are almost always considered. The first is the nonlinearity at the heart 
of the Jaynes-Cummings model, arising when each cavity embeds a two-level system, and 
gives rise to the model known as the Jaynes-Cummings- Hubbard (JCH) model [7]. The 
second is the standard third order optical nonlinearity. It can arise from the nonlinear 
optical response of the cavity material [8], or from polariton-polariton interactions in 
the case of exciton-polaritons [9] or polaritons originating from driven few-level atoms 
[TU1 [TTT IT2l IT3"] . In this second case, the resulting Hamiltonian is sometimes called the 
Kerr-Hubbard (KH) Hamiltonian, as it corresponds to the Bose-Hubbard Hamiltonian 
with the addition of driving fields and dissipation. 

A vast theoretical literature has appeared in the last decade, describing several 
ways in which quantum phase transitions could be observed in coupled cavity arrays 
(as reviews, see Refs. pHl EH])- Early works [151 EH1 El EH EE] have focused on 
demonstrating how specific optical systems could be mapped onto the JCH or KH 
Hamiltonians, thus displaying typical critical phenomena of strongly correlated systems 
such as the emergence of a Mott insulator or Tonks-Girardeau phase. These works 
also suggested that the two models can coincide under specific conditions [TTIEH]- The 
case of disordered arrays was also considered [20J, where the quantum phase transition 
between a superfluid and a Bose-glass phase was characterized. 

In these works however the driven-dissipative nature of the optical systems under 
investigation was only marginally taken into account, while an assumption of a lossless 
system at equilibrium was essentially employed to model the occurrence of quantum 
phase transitions. Within this assumption, the typical condition for the onset of strong 
quantum correlations is, analogously to the prototypical case of the Mott insulator, that 
U ^> J, where U represents the on-site interaction energy per particle and J the hopping 
energy between adjacent sites. It was only in subsequent works [91 [211 [221 [23] that the 
important role of dissipation was clarified, showing that driven-dissipative systems may 
differ considerably from their lossless counterparts. 
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In a driven-dissipative system, each photon is confined in a cavity only during a 
time lapse corresponding to the cavity lifetime, after which it eventually radiates out of 
the cavity array. An external source - typically one or more laser beams - replenishes 
the system with new photons, and both time-dependent and steady-state regimes can be 
achieved depending on the experimental protocol. General wisdom suggests that, in a 
driven-dissipative system, quantum correlated states of many particles can be obtained 
only if the particle lifetime is long enough. More precisely, if one particle is added to 
the system, the time it takes to acquire quantum correlations with the other particles 
present in the array must be much shorter than the time the particle spends inside 
the system. The rates of these two processes are characterized by the effective on-site 
interaction energy per particle U, and by the loss rate 7 respectively. It is clear then, 
that the condition U/j 1 is required in order for the many-particle state to acquire 
quantum correlations. The most illustrative example of this concept can be found 
already in the analysis of the single site KH problem, as was done for example in the 
context of polaritons in micropillars [9] , where it is shown that the condition U/7 > 1 
is required in order to reach the blockade regime [121 Q2] an d emit photons with sub- 
Poissonian statistics. The parameter {7/7 is an example of cooperativity [21] and it is 
indeed related to the atomic cooperativity parameter that enters the dissipative Jaynes- 
Cummings model of CQED [17] (although with a caveat [19]) or more specific polaritonic 
models of few-level atoms pi]]. According to this intuitive argument, systems with small 
cooperativity should be restricted into non-entangled quantum states, generally having 
a classical analog. 

Reaching high cooperativity in realistic systems is very challenging. Currently, the 
state-of-the-art value of this parameter is of the order of 1 in both atom-semiconductor- 
resonator [25] and all-semiconductor systems [26], while it reaches values of about 10 in 
atom-cavity systems [27] and of about 50 in superconducting circuits [28]. 

The high cooperativity argument however relies on a rather strong assumption, 
namely that at the initial time an extra particle added to the system is set into a 
separable state, namely is not already entangled to the other particles. While this 
would be the case for atomic gases, optical systems have as a unique feature that the 
source is a coherent field, with controlled amplitude and phase that are passed on to 
the many-photon system inside the cavity array. Quantum interference effects might 
then produce quantum correlations directly as a result of the driving process, in a time 
shorter than the lifetime. Following this general idea, we have recently shown that 
quantum effects such as sub-Poissonian statistics [29J can be expected in a system with 
very low cooperativity, under particular conditions. This is possible due to quantum 
interference between different excitation pathways and its sensitivity to small energy 
shifts caused by the interaction of pairs of particles [30]. Theoretically, the quantum 
interference mechanism can also be exploited to create bipartite entangled states in a 
system of three coupled cavities [31], where the entanglement is based on the amplitude 
and phases of quantum optical fields. This continuous variable (CV) entanglement 
[32] can be detected using a homodyne interference setup [33], not being conditioned 
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on single photon detection. Being issued of a quantum interference process, differently 
from the blockade effect, this mechanism is rather sensitive to the rate of pure dephasing 
characterising each quantum resonator. It is then very important to quantitatively 
account for pure dephasing in all model calculations. 

A remaining question is whether multipartite entangled states can be produced in 
a larger array of linearly coupled optical resonators in the regime of low cooperativity. 
Besides the fundamental interest, multipartite entanglement represents an important 
resource for applications in quantum information science [SI EE] . While it is possible 
to produce multipartite states by coupling parametric down-conversion processes in 
nonlinear crystals [361 E3 EEl EHJ HQ] or using optical frequency combs [UJ 02], arrays 
of optical cavities can be realized as compact solid-state systems, which are more 
appropriate for devices. 

Here, after reviewing our analysis of a three-cavity array for the generation of two 
entangled photons [31], and extending it to the JCH Hamiltonian, we study middle- 
sized arrays of up to eight coupled modes. We demonstrate two schemes for the gen- 
eration of multipartite CV-entanglement, which is assessed via the violation of van 
Loock-Furusawa inequalities [13] (a multipartite CV analogue of Bell inequalities). Each 
scheme is characterized by an array topology and by the selection of the subset of cavities 
that are not laser-driven in a mutually coherent fashion. Four-party CV-entanglement is 
achieved at an optimal amplitude of the driving field, corresponding to an average mode 
occupation of slightly below 1, and is tested against a realistic value of the pure dephas- 
ing rate. By modifying the relative phases of the driving fields, it is shown that some 
control of the entanglement is possible, holding promise for the realization of devices 
for the controlled generation of multipartite entangled photons. Previous theoretical 
calculations of small arrays have relied on the direct numerical solution of the Von Neu- 
mann equation for the density matrix in a truncated Fock space (particle number basis). 
While this is an accurate approach within quantum optics, it is numerically heavy and 
cannot be applied to systems of more than a few coupled cavities. In this work, we make 
use of the wavefunction Monte Carlo (WFMC) technique [S], which gives access to the 
exact solution of the Von Neumann equation by trading off speed against low memory 
consumption. 



2. Theory 

We consider arrays of coupled nonlinear optical resonators under near-resonant 
excitation. The theoretical model will be suitable for the description of a variety 
of systems that can be divided into two types. First, polaritonic systems, in which 
the polaritons are well described as bosons with Bose-Hubbard-like interactions. This 
includes not only exciton-polaritons in semiconductor microcavities, but also the 
polaritons that arise whenever many two-level systems are resonantly coupled to an 
optical resonator [TO], [121 CBJ- Second, the case of a single few- level system (e.g., an 
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atom or quantum dot) coupled to an optical resonator. In this latter case, the Kerr 
nonlinearity is a good description in the limiting case of the dispersive regime [TTf [T9] . 

Under near-resonant excitation, higher energy modes in realistic systems can be 
neglected, such that in the model each resonator is described by a single mode with 
energy Ej. The cavities are also defined by a constant U characterizing the nonlinear 
energy per particle, which we assume equal in each cavity for simplicity. A matrix 
Jjk describes the tunnelling rates between pairs of resonators, which are related to the 
topology of the system (i.e., the spatial arrangement of the cavities). Under these 
assumptions the KH Hamiltonian [15] of the system is: 



where hj are the Bose annihilation operators associated with the modes and Fj are 
the complex amplitudes of driving fields originating from the optical pumps acting 
on the modes. Note that this Hamiltonian is written directly in the rotating frame 
of the driving field, so that Fj is a constant in time and the energies Ej are expressed 
relative to the optical pump energy fvjjQ. For simplicity, we are assuming monochromatic 
excitation, namely that the optical fields acting on each cavity have the same energy. If 
this condition is lifted, then typically the system does not admit a steady state solution 
even for stationary driving fields, The study of this situation lies beyond the scope of 
the present work. 

The description of a driven-dissipative system requires the introduction of a 
coupling to an external reservoir. Here we assume two effects originating from this 
coupling. First, to each cavity mode a dissipation process Out ci r£ite 7 is associated, 
corresponding to the finite lifetime of a photon inside the cavity. Second, we also 
introduce a pure dephasing rate associated to each mode, at a rate T. Pure dephasing 
is the result of an interaction of the photon with the random environment of the 
bath without leakage of the photon out of the cavity (for example the exciton-phonon 
scattering mechanism in the case of exciton-polaritons in semiconductor microcavities 
[32] )■ In terms of the density matrix of the system, the pure dephasing describes 
a decay of the off-diagonal terms of the density matrix without a corresponding 
decay of the diagonal ones. For simplicity, all modes are assumed physically similar, 
being characterized by common dissipation and dephasing rates. Dissipation and pure 
dephasing can be accounted for in the Von Neumann equation for the density matrix by 
introducing corresponding Lindblad terms [36] , as detailed in the next section. We note 
that the weakly nonlinear regime considered in this work is characterized by 7 ^> U and 



We will consider two different approaches for studying the dynamics of this driven- 
dissipative system: the direct numerical solution of the Von Neumann equation for the 
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density matrix, and the Wavefunction Monte Carlo (WFMC) approach. The former is 
accepted as the most accurate within the framework of quantum optics, however is not 
feasible for large systems due to the exceedingly large size of the corresponding density 
matrix. The latter has been shown to give equivalent results and is suitable for larger 
systems [H]. We will provide numerical evidence of the equivalence between the two 
approaches when reviewing our recent result in the case of a system of three coupled 
modes [31] where bipartite entanglement is expected. Afterwards, we will consider 
systems of seven and eight coupled cavities and demonstrate that they can give rise 
to multipartite entanglement of four photons under appropriate conditions. For these 
larger systems, we will restrict to the use of the WFMC approach, as a direct solution 
of the density matrix equations would be numerically unfeasible. 



2.1. Von Neumann Equation 

The quantum optical behaviour of our system can be fully described using the Von 
Neumann equation for the density matrix, p: 



+ i- ^2 \ 2 a jpOj - a)ajp - pa^Ojj 

j 

+ ( 2 %^ - ™?P - P^j) > ( 2 ) 

3 

The second and third terms on the right-hand side of this equation are Lindblad 
terms required to describe the Markovian coupling between the system and a random 
environment. They describe respectively the dissipation and pure dephasing processes 
|46j . Here, nj = d^dj is the occupation operator of the j-th mode. Equation [2] can be 
solved by expanding the density matrix over a particle number basis in a similar way to 
that done in Ref. [9]: 

p = |m, n 2 , n N ) (n[,n' 2 , ...,n' J v|Pn 1 ,n 2 ,...,n JV X 1 ^,...,n' lV (3) 

It is necessary to truncate at a given particle number, m < n max , which is accurate in 
the regime of small occupations. Substitution of Eq. [3] into Eq. [2] gives a finite set of 
evolution equations for the elements of the density matrix, p ni ,n 2 ,...,njv ! rei,n! ! ,...,n5 v ( see the 
supplemental material of Ref. [29] for more detail). One can then time evolve the system 
numerically from the vacuum to obtain the steady state density matrix, from which all 
observable quantities can be calculated. It is often easier however to solve directly for 
the steady state - when assuming constant Fj's - by solving the linear homogeneous 
set of equations obtained for dp/dt = 0, with the additional condition tr(p) = 1. 
As an additional simplification of the numerical calculations, it is straightforward to 
further truncate the basis of the wavefunction such that only states \n\,n2, nj\r) with 
J2 m n m. < n max are included. For the same value of n max , this prescription amounts to 
selecting a considerably smaller Hilbert space. The justification of this choice lies in the 
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fact that this latter truncation scheme retains all states within a maximum total energy, 
responding to an intuitive perturbation criterion, while the former one includes many 
more states all characterized by higher energy, which are presumably less relevant in a 
perturbation expansion. 



2.2. Wavefunction Monte Carlo Technique 

The WFMC technique was developed and reported in Ref. [H]. Here we give a brief 
summary of the basic elements applied to the class of systems we are investigating. 
Time evolution of the density matrix using the Von Neumann equation is replaced with 
time evolution of a wavefunction, \ip(t)), using the non-Hermitian Hamiltonian 

ih v— v „ ih„\-^ „ 
H = H s --l^n--T^l (4) 

n n 

where n n = a) n a n . Since this Hamiltonian is non-Hermitian it is important to 
renormalize the wavefunction at every timestep in a numerical calculation such that 
the normalization condition (ip(t)\ijj(t)) = 1 is preserved during the time evolution. 

In addition to evolution with the non-Hermitian Hamiltonian, there is a possibility 
of a stochastic quantum jump in the wavefunction between the times t and t + St. 
Two different types of jumps are possible, associated with either the dissipation or pure 
dephasing: 

\^(t + 5t)) =a m \^{t)), (5) 
W + 5t)) =n m \^(t)), (6) 

which occur, respectively with probabilities 

Sp m (t) =7<ft(V(*) \n m \^(t)) (7) 

6p' m (t) = T5t(m\n 2 m \m)- (8) 

Using these probabilities, the WFMC procedure is equivalent to the simulation of a 
stochastic process. By repeating the calculation over a large enough ensemble of time 
sequences of stochastic quantum jumps, a distribution of different wavefunctions can 
be obtained at each time. Note that one should choose a short enough timestep such 
that Yl m (dp™ + 5Pm) ^ 1- It is also important to renormalize the wavefunction after 
a jump, since Eqs. [MS do not preserve the normalization. 

In Ref. [H] it was shown that the WFMC procedure is equivalent to the Von 
Neumann equation approach when observable quantities are calculated by averaging 
over multiple wavefunction realizations. The expectation value of the operator O is 
given by averaging over Nr different realizations of the wavefunction (t)^, each 
obtained from a different time sequence of quantum jumps: 

N R 

[0)^ w ^(^(t)\0\^(t)) (9) 
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In this paper, the systems that we study admit steady state solutions. After a steady- 
state is reached one can also perform time averaging, which reduces the statistical error 
of the technique. 

The advantage of the WFMC technique is that the wavefunction is expanded over a 
basis with less than (n max + 1) M states - where M is the number of cavities - in contrast 
to the direct solution of the equation for the density matrix where (n max + l) 2M equations 
need to be solved. This enables the solution of problems that would otherwise require 
amounts of memory too large to be handled with standard modern computers. As for 
the direct solution of the master equations, one can reduce the size of the basis by 
including only states with ^2 m n m < n max . In all numerical solutions of the problem 
presented below, we have carefully checked the convergence of the result as a function 
of n max - Typically, convergence was reached for values of n max ranging from 8 to 12 
depending on the array and parameter range under study. 

3. Two Party Entanglement in a Three Coupled Mode System 

In this section we review our recent result [31] on a system of three cavities coupled 
in a chain, where the central cavity is driven by a near-resonant optical excitation (see 
inset in Fig. [T^i). This system is sufficiently simple to be treated using both the full 
Von Neumann equation approach (section l2.ip and the WFMC approach (section l2.2p . 
which allows to test the effectiveness of the WFMC method. In the system under 
consideration, it has been previously shown that two-party CV entanglement [32] can 
develop between the two non-driven modes due to the presence of quantum interference 
in the system [3T] . 

In analogy to Bell's result for discrete variable entanglement, CV entanglement 
between modes n and m is characterized by the violation of an inequality pTH HE] : 

1 < S nm = V (p n - p m ) + V (q n + q m ) , (10) 

where the amplitude and phase operators are defined as p n = [a n + a n )/2 and 
Qn — { a n — a V) /(2i), respectively. The variance of an operator is defined by V{0) = 
(O 2 ) — (O) 2 . To be sure that inequality [TUl evidences entanglement if it is present, it 
may be necessary to minimize the value of S13 over the phase references of the different 
modes. Local operations allow to make the transformation a n 1— > a n e~ lr ^ n and the 
minimum value of S nm is given by: 

S nm = 1 + (a n a n ) + (a^cim) - (a n )(a n ) - (a,| n )(a m ) 

-2\J (aWm) - (at) (aL) y/ {a n a m ) - (a n )(a m ) (11) 

The expectation values entering this equation can be calculated from the density matrix, 
p, or from Eq. |9l depending on the theoretical approach used to model the system. 

The variation of §13 with driving field amplitude is shown in Fig. (TJ using the 
same set of parameters as in Ref. [31]. These parameters were originally chosen to 
correspond to exciton-polariton boxes [5] of size 3/im separated by 1/zm, for which 
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Figure 1. (color online) a) Variation of the parameter S13 with driving field amplitude 
for a chain of three coupled cavities. Spots correspond to results from the Von Neumann 
equation approach, while curves correspond to results from the WFMC approach. 
Different traces (spots and curves) correspond to different values of pure dephasing, T 
(values marked on plot). The inset shows the topology of the cavity array, where the 
central cavity is driven, b) % error of WFMC approach, c) Variation of the occupation 
numbers of the cavities with driving field amplitude. Again, spots correspond to results 
from the Von Neumann equation approach, while curves correspond to results from 
the WFMC approach. The large spots correspond to the average occupation of the 
modes (ni) = (713), while the small spots correspond to (712) - 



J = 0.5meV [29] can be obtained from solution of the Schrodinger equation in a two well 

potential (this value is also in agreement with experimental measurements [4~9|l50]). It is 
well-known how to calculate the nonlinear interaction strength [9] and we take the value 
U = 0.012meV in agreement with recent experimental data [5TJ [52j [53]. The exciton- 
polariton decay rate was taken as 7 = 0.044meV, which was measured in Ref. [M] and 
we chose E1 + E3 = — 0.06meV (Ei = E 3 in a and b), E 2 = 0.08meV. The slight detuning 
between the modes E\$ and E2 was previously found to give the smallest value of S13 
for fixed J and 7 in Ref. [31]. It should be noted that the master equations issued of 
the Hamiltonian (CQ), with dissipation and dephasing rates 7 and T respectively, can be 
reduced in a dimensionless form by measuring all energy parameters with respect to 7. 
Hence, all results shown here are invariant under rescaling of all energies and therefore 
hold for systems characterized by different ranges of parameters such as, for example, 
superconducting devices [28] . 

In agreement with the calculations based on the Von Neumann equations (spots 
in Fig. [1]), the WFMC approach (curves) demonstrates entanglement between the 
first and third modes in the system via violation of the inequality 1 < S13. The 
differences between the two approaches are shown in Fig. [lb, which illustrates that 
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the WFMC approach is an accurate alternative technique for our problem. It is 
notable that the technique becomes slightly less accurate for large values of dephasing, 
where entanglement is lost. However, since realistic estimates of the dephasing rate in 
semiconductor micropillars are in the range of tenths of /ie V [H] , the technique remains 
very accurate for the purposes of this paper. A minimum in the value of S13 is obtained 
for a driving rate of F ~ 0.7 meV, while for larger rates this quantity rapidly increases 
to values above 1, thus recovering a non-entangled state as expected in the limit of high 
average occupation of the modes. The variation of the occupation numbers with driving 
field amplitude is shown in Fig. [lb, where again there is full agreement between the two 
approaches. It is important to notice that entanglement occurs for average occupation 
below unity, where nonclassical effects are indeed expected. 




F (meV) F (meV) 

Figure 2. (color online) a) Variation of the parameter S13 with driving field amplitude 
for a chain of three coupled cavities, in the case of the JCH model. Different traces 
correspond to different values of pure dephasing acting on the two-level systems, T s 
(values marked on plot), b) Variation of the occupation numbers of cavities 1 and 3 
(full lines) and of cavity 2 (dashed lines) with driving field amplitude. Different colors 
denote different values of T s as in panel a). 

To conclude this section on bipartite entanglement, we extend our previous result 
[3Tj by modelling the same array geometry in which the Kerr-Hubbard Hamiltonian 
is replaced by a Jaynes-Cummings-Hubbard one. The question whether the present 
mechanism holds in the JCH case is a nontrivial one, as the JCH Hamiltonian is not 
in general equivalent to a KH one [T7] [T9] . In the limit of very large detuning between 
the cavity and the two-level system, the JCH case can indeed be associated to an 
effective KH Hamiltonian. The case of major relevance however is the one close to 
resonance, as in all atom-cavity, semiconductor-QD, or Josephson-qubit systems. The 
weak cooperativity limit was already studied in the case of two coupled cavities, one of 
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which coupled to a two-level system [20], showing that an effect quantitatively analogous 
to the one arising in the KH case can be obtained for moderate cooperativity (g/'y — 1, 
where g is the Jaynes-Cummings coupling constant introduced below). Here we study 
the three-cavity configuration, assuming that each of the two side ones is coupled to a 
two-level system. The Hamiltonian of this model system is 



where &j± are the raising and lowering operators of the j-th two-level system respectively. 
Lindblad terms for the two-level systems are written analogously to those for the cavity 
modes. Here, in addition to taking a single value for J and for g, we also assume 
degenerate side cavities, namely E\ = E% and E\ s = E 3s . For the numerical calculations, 
we take a dissipation rate of the two-level systems 7 S = 7/IOO. This value should 
model nonradiative decay channels as well as decay into photonic modes other than 
the cavity one, and matches typical values extrapolated for semiconductor systems 
[55] . We neglect pure dephasing processes for the three cavity modes and retain only a 
pure dephasing rate T s for the two-level systems, as this process is expected to be the 
dominant source of dephasing. Fig. [2^ shows the variation of the parameter S 13 with 
driving field amplitude, computed as in the KH case, while in Fig. [2b the corresponding 
average occupations are plotted. Here we assumed parameters as in the system of a 
quantum dot in a photonic crystal cavity studied in Ref. [26], namely g = 141 y^eV and 
7 = 53 /ieV. We have run a minimum search procedure on the energies and found an 
optimal violation of the generalized Bell's inequality for values E\ — E$ — 0.06 meV, 
E 2 = -0.12 meV, E Xs = E 3s = 0.22 meV, and J = 0.29 meV, as shown in Fig. [2k, At 
the optimal value of F the violation is of the same order as in the KH case, showing 
that the mechanism studied in this work does not necessarily rely on a specific kind of 
nonlinearity. Curves computed for the same values of the parameters but in presence 
of a finite T s are also shown. The pure dephasing of an InGaAs quantum dot has 
been extensively characterized by coherent four-wave mixing measurements, and its rate 
typically amounts to T s = 1 /ieV [56], or even smaller [57]. For this value we see that 
the value of S13 is only slightly affected and the non-separability of the steady state is 
robust to pure dephasing rates up to ten times larger. Coupling between photonic crystal 
cavities is currently feasible and can be easily tuned [58] . This result thus suggests that 
current semiconductor technology would be advanced enough to produce a proof-of- 
principle device for bipartite continuous-variable entangled photon generation. 



3 




(12) 
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4. Multiparty Entanglement 

The scheme discussed in the previous section can be generalized to many resonators. 
The objective is the generation of cluster-entangled states, that are a universal resource 
for quantum information processing [34]. In very simple terms, a cluster state is a 
state of many qubits arranged on an array, such that each qubit is entangled with the 
neighbouring ones. It is then natural to expect the emergence of a cluster state from the 
spatial repetition of the three cavity system discussed above. Hence, as a first example 
of an array topology able of generating such a multipartite entangled state, we consider 
eight cavities arranged on a ring, with next neighbour coupling, as illustrated in the 
inset to Fig. [3k. Alternate cavities are optically driven, so as to repeat the three-cavity 
pattern of the previous section. Quadripartite entanglement should build up between 
the four non-driven modes. For simplicity we assume that the driving fields have equal 
amplitudes, such that: -F evcn = F and F odd = 0. As before, we compute the steady state 
solution of the equations for the density matrix, issued from the Hamiltonian (CQ), in 
presence of dissipation and pure dephasing. 




F (meV) F (meV) 

Figure 3. (color online) a) Variation of the quadripartite entanglement parameter / 
with driving field amplitude in the case of a ring of eight coupled modes with near- 
resonant optical pumping of alternate modes (illustrated by the inset). The different 
curves correspond to different values of pure dephasing (marked on the plot), b) 
Variation of the average occupation of the modes (note that all modes have the same 
average occupation due to symmetry). Parameters: T = 0.044meV, U = 0.012mcV, 
J = 0.16meV. The mode detunings with respect to the pump energy were chosen to 
minimize /: E oc id = — 0.03meV, E even = O.llmeV. 

Quadripartite continuous variable entanglement can be evidenced by the 
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simultaneous violation of three inequalities [ID]: 

h = V(p - pi) + V(g + qi + 92^2 + gm) > 1, (13) 
h = V(pi - P2) + V(g qo + gi + g 2 + ^3^3) > 1, (14) 
h = V(p 2 - pa) + V(#o<?o + + <?2 + $3) > 1, , (15) 
where are arbitrary real parameters that should be chosen to minimize the largest 
value of the set of I n , that is, I = min 9i (max{Ii, I2, ^3})- As in the bipartite case, an 
additional optimization is made by varying the phase reference of each of the modes. 

Results from the WFMC approach are shown in Fig. [3h, where the variation of 
/ with driving field amplitude is shown for a range of values of pure dephasing. For 
small dephasing, the value of I drops below 1 indicating the presence of quadripartite 
entanglement between the modes Si, 63, 65, and 07. This entanglement survives realistic 
values of pure dephasing in the range of tenths of /xeV [J5]. The entanglement is lost 
for large driving field amplitude, as before, because the system approaches a classical 
behaviour for large average occupation of the modes. For the optimum value of /, the 
average occupations of the modes, shown in Fig. [3b, are approximately 0.45 each, again 
well in the quantum regime. 




F (meV) F (meV) 

Figure 4. (color online) a) Same as in Fig. [3] for a linear chain of modes (illustrated 
by the inset in a). The mode detunings with respect to the pump energy were chosen 
to minimize /: E Q dd = — 0.02meV, E cven = — 0.19meV. The other parameters were the 
same as in Fig. [3] In (b) the solid curves show the values of (rii) = (n?), while the 
dashed curves show the values of (723) = (n^) (these mode occupations are equal due 
to symmetry). 

Figure H] shows an alternative scheme of operation based on a chain of seven modes 
(illustrated in the inset to Fig. Hk). Similar to the case in Fig. [31 driving alternate modes 
allows quadripartite entanglement of the non-driven modes. This demonstrates the 
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potential of arrays of weakly nonlinear cavities to generate entanglement in a variety of 
topologies that can be physically engineered. Again, the optimum value of / corresponds 
to an average mode occupation of approximately 0.45 (see Fig. Hb). 

It is possible to switch off the quadripartite entanglement by changing the 
configuration of the driving fields. Removal of one of the driving fields in the eight 
mode system (F 8 = 0) causes the loss of quadripartite entanglement as shown by the 
dashed curve in Fig. [5^. Note that this system, with one driving field removed, is 
not identical to the chain of seven modes considered in Fig. HI With one driving field 
removed, and no re-optimization of parameters, the occupations of the different modes 
in the system may be very different, as shown in Fig. [5b. Indeed, when the mode 
occupations are not similar, one does not expect a strong entanglement since there is a 
dominant probability of detecting particles in a subset of modes. 




F (meV) F (meV) 

Figure 5. (color online) a) The effect of removing one of the driving fields in the 
eight mode ring scheme is shown by the dashed curve (the solid curve is the same 
as in Fig. [3|)a. b) The variation of the average mode occupations with all eight 
driving fields (solid curve) and one driving field removed (dashed curves). Due to 
symmetry, (71,3} = (77,5) and (ni) = (n?) when the driving field F% is removed. The 
other parameters were the same as in Fig. [3^,, with T = 0. 

In addition, the entanglement is sensitive to the relative phase of the driving fields. 
Keeping all fields switched on but varying the phase of one of the fields (F 8 ) is also able 
to switch off the entanglement, as shown in Fig. [6^ where the entanglement parameter is 
plotted for a range of relative phases. Again, there is an associated imbalance of average 
mode occupations when varying the phase of one of the fields as shown in Fig. [6b. 

It is also interesting to study the behaviour of the bipartite entanglement parameter, 
S nm , between two selected modes among the four being studied. Figure [7^ shows the 
variation of S nm , evaluated between pairs of signal modes, for the eight mode ring with 



Multimode entanglement in coupled cavity arrays 



15 




0.90 



0.0 



0.1 0.2 0.3 

F (meV) 







0.1 



0.2 



0.3 







F (meV) 



Figure 6. (color online) The effect of varying the phase of one of the driving fields 
(numbers marked on the plot indicate the relative phases used) . In (b) the solid curves 
show the values of {h\) = {fir), while the dashed curves show the values of (773) = (7x5) 
(these mode occupations are equal due to symmetry). Parameters the same as in 
Fig. EJa,, with T = 0. 



all driving fields switched on (the case of Fig. [3k,). It can be seen that there is no 
bipartite entanglement between any pair of modes since inequality [10] is not violated - 
despite the fact that quadripartite entanglement exists, this entanglement is locked in 
the four entangled modes. This is similar to the property of the Greenberger-Horne- 
Zeilinger (GHZ) three-qubit state: tracing out one party leaves the system unentangled. 
Figure [7b shows the variation of S nm for the case when one driving field is removed from 
the eight mode ring. In this case the system displays neither quadripartite nor bipartite 
entanglement. 

5. Conclusion 

Solid-state systems of arrays of weakly nonlinear coupled cavities represent a viable 
route toward the generation of quantum correlated phenomena in a compact device. 
In this work, we have demonstrated that the emergence of quantum correlations is 
not necessarily restricted to systems with high cooperativity - namely large nonlinear 
energy as compared to the dissipation rate - but can occur in the weakly nonlinear 
limit thanks to quantum interference between different pathways of coherent excitation. 
Using the Wavefunction Monte Carlo approach we were able to study the quantum 
optical behaviour of moderately large systems. In particular, we have reviewed our 
early results for the system of three coupled cavities with a driving field acting on 
the middle one [31]. This system may be viewed as the elementary building block 
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Figure 7. (color online) a) Variation of the bipartite parameter S13 with driving field 
amplitude for the eight mode ring for different values of pure dephasing. Note that 
S13 ~ Sis and other combinations of modes are equivalent to either S13 or 5x5 due 
to symmetry, b) Variation of the bipartite parameters in the eight mode system with 

^8 = 0. 

for the design of more complex arrays for multipartite entanglement generation. We 
have also extended our previous study by modelling the same array in which the Kerr 
nonlinearity is replaced by the Jaynes-Cummings one. Also in this case, we have 
obtained a significant violation of the generalized Bell's inequality, suggesting that this 
entanglement generation paradigm can be ported to CQED systems such as photonic 
crystal cavities embedding semiconductor quantum dots [26]. As an example of a 
system able of generating multipartite entanglement, we proposed an array of eight 
coupled cavities, which could be realized in a variety of structures including photonic 
crystals, atom-cavity systems and semiconductor micropillars. We have shown that 
the entanglement is dependent on the configuration of driving fields and that it can 
also be generated in systems with other topologies, such as a linear chain of seven 
modes. This result represents a very important proof of principle for the control of the 
entanglement, where one would be able to alter the system topology by electric [52] 
or magnetic [60] fields as well as vary the pattern of the incident optical field. The 
system has more potential for scalability and control than previous schemes based on 
parametric scattering [61]. 

Future work should have, as an objective, to find a deterministic link between 
the type of entangled state on one side, and the array topology and optimal value 
of parameters on the other side, similarly to what has been done for the elementary 
system of two coupled cavities [30] . The full numerical solution of the steady-state Von 
Neumann equation for the density matrix is a very demanding computational task and 
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the arrays studied here are the largest systems that can be tackled by exact numerical 
methods. The study of larger systems will require the adoption of approximated methods 
able of correctly describing the quantum correlated nature of the states. Matrix product 
operators could be the election method for one- dimensional arrays [62| [63] , while arrays 
of higher dimensions would require a generalized tensor network approach [M] • Finally, 
further investigations are still required to assess the usefulness of the entangled states 
within the field of quantum technology. Indeed, schemes of quantum computation 
based on multipartite entangled states [Hj remain promising and have been adapted for 
continuous variable systems [35]. However, the criteria for multipartite entangled states 
to be useful resources have not been completely developed. 
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